# New descriptive visualisation for SCST paper ---------------------------------

# Set up survey design for srvyr package ---------------------------------------
svypool <- pool %>%
  as_survey_design(id=hh, strata=village, weight =weights, nest = TRUE)

# Access means grouped by state ----------------------------------------------

access_means <- svypool %>% 
  group_by(state, wave, scst) %>% 
  summarize("Electricity (Grid)" = survey_mean(usegrid, na.rm = T),
            "LPG" = survey_mean(uselpg, na.rm = T)) %>% 
  mutate(wave = factor(wave, labels = c("0" = "2014","1" = "2018")),
         scst = factor(scst, labels = c("0" = "All Others", "1" = "SC/ST")),
         state = case_when(state == 9 ~ "Uttar Pradesh",
                           state == 10 ~ "Bihar",
                           state == 19 ~ "West Bengal",
                           state == 20 ~ "Jharkhand",
                           state == 21 ~ "Odisha",
                           state == 23 ~ "Madhya Pradesh")
         ) %>% 
  gather(attribute, mean, -c(state, wave, scst)) %>% 
  tbl_df()

access_means %>% 
  filter(!grepl(attribute, pattern = "se")) %>% 
  ggplot(aes(x = wave, y = mean, colour = scst, group = scst)) + 
  geom_line() +
  geom_point() +
  facet_grid(attribute~state) +
  labs(y = "Proportion of households with access",
       x = NULL,
       colour = "Social Group") +
  scale_y_continuous(labels = scales::percent) +
  theme_bw() +
  theme(strip.background = element_rect(fill = NA, color = "black"),
        legend.position = "top")

ggsave(here("Manuscript","Figures","descr_acc.png"),
       device = "png", width = 8, height = 6*0.618, units = "in")

# Electricity supply duration means grouped by state -------------------------

elechrs_means <- svypool %>% 
  filter(usegrid == 1) %>% 
  group_by(state, wave, scst) %>% 
  summarise(day = survey_mean(duration_day, na.rm = T),
            eve = survey_mean(duration_night, na.rm = T)) %>% 
  mutate(wave = factor(wave, labels = c("0" = "2014","1" = "2018")),
         scst = factor(scst, labels = c("0" = "All Others", "1" = "SC/ST")),
         state = case_when(state == 9 ~ "Uttar Pradesh",
                           state == 10 ~ "Bihar",
                           state == 19 ~ "West Bengal",
                           state == 20 ~ "Jharkhand",
                           state == 21 ~ "Odisha",
                           state == 23 ~ "Madhya Pradesh")) %>% 
  tbl_df()

day <- elechrs_means %>% 
  ggplot(aes(x = wave, y = day, colour = scst, group = scst)) + 
  geom_line() +
  geom_point() +
  facet_wrap(~state) +
  labs(title = "Daily supply hours (0-24 hours)",
       y = "State mean supply hours",
       x = NULL,
       colour = "Social Group") +
  theme_bw() +
  scale_y_continuous(breaks = seq(0,24,4),
                     limits = c(0,24)) +
  theme(strip.background = element_rect(fill = NA, color = "black"))

eve <- elechrs_means %>% 
  ggplot(aes(x = wave, y = eve, colour = scst, group = scst)) + 
  geom_line() +
  geom_point() +
  facet_wrap(~state) +
  labs(title = "Evening supply hours (0-6 hours)",
       y = "State mean supply hours",
       x = NULL,
       colour = "Social Group") +
  theme_bw() +
  scale_y_continuous(breaks = seq(0,6,2),
                     limits = c(0,6)) +
  theme(strip.background = element_rect(fill = NA, color = "black"))

# Electricity supply disruption means grouped by state ------------------------

elecdays_means <- svypool %>% 
  filter(usegrid == 1) %>% 
  group_by(state, wave, scst) %>% 
  summarise(outages = survey_mean(reliability, na.rm = T),
            lowvolt = survey_mean(quality_low, na.rm = T),
            fluct = survey_mean(quality_fluc, na.rm = T)) %>% 
  mutate(wave = factor(wave, labels = c("0" = "2014","1" = "2018")),
         scst = factor(scst, labels = c("0" = "All Others", "1" = "SC/ST")),
         state = case_when(state == 9 ~ "Uttar Pradesh",
                           state == 10 ~ "Bihar",
                           state == 19 ~ "West Bengal",
                           state == 20 ~ "Jharkhand",
                           state == 21 ~ "Odisha",
                           state == 23 ~ "Madhya Pradesh")) %>% 
  tbl_df()

outages <- elecdays_means %>% 
  ggplot(aes(x = wave, y = outages, colour = scst, group = scst)) + 
  geom_line() +
  geom_point() +
  facet_wrap(~state) +
  scale_y_continuous(limits = c(0,6)) +
  facet_wrap(~state) +
  labs(title = "Monthly outages",
       y = "State mean outage days",
       x = NULL,
       colour = "Social Group") +
  theme_bw() +
  theme(strip.background = element_rect(fill = NA, color = "black"))

lowvolt <- elecdays_means %>% 
  ggplot(aes(x = wave, y = lowvolt, colour = scst, group = scst)) + 
  geom_line() +
  geom_point() +
  facet_wrap(~state) +
  scale_y_continuous(limits = c(0,6)) +
  labs(title = "Monthly Low Voltage Days",
       y = "State mean low-voltage days",
       x = NULL,
       colour = "Social Group") +
  theme_bw() +
  theme(strip.background = element_rect(fill = NA, color = "black"))

fluct <- elecdays_means %>% 
  ggplot(aes(x = wave, y = fluct, colour = scst, group = scst)) + 
  geom_line() +
  geom_point() +
  facet_wrap(~state) +
  scale_y_continuous(limits = c(0,6)) +
  labs(title = "Monthly Voltage Fluctuation Days",
       y = "State mean voltage fluctuation days",
       x = NULL,
       colour = "Social Group") +
  theme_bw() +
  theme(strip.background = element_rect(fill = NA, color = "black"))

day / (eve + outages) / (lowvolt + fluct) + plot_layout(guides = "collect") & theme(legend.position = 'bottom')

ggsave(here("Manuscript","Figures","descr_elec_all.png"),
       device = "png", width = 8, height = 12, units = "in")

# Energy cost burden visualisation ---------------------------------------------

energy_costs <- pool %>% 
  transmute(hh = hh,
            state = state,
            wave = ifelse(wave == 0, "2014", "2018"),
            scst = ifelse(scst == 1, "SC/ST", "All Others"),
            month_exp = month_exp,
            fuel_cost = fuel_cost,
            elec_cost = rowSums(select(.,grid_spend, mg_spend, kero_spend), 
                                na.rm = T)) %>% 
  mutate(elec_share = round(elec_cost / month_exp,3),
         fuel_share = round(fuel_cost / month_exp,3),
         total_share = round((elec_cost + fuel_cost) / month_exp,3),
         state = case_when(
           state == 9 ~ "Uttar Pradesh",
           state == 10 ~ "Bihar",
           state == 19 ~ "West Bengal",
           state == 20 ~ "Jharkhand",
           state == 21 ~ "Odisha",
           state == 23 ~ "Madhya Pradesh"
         )) %>% 
  filter(month_exp > 0)

energy_costs %>% 
  group_by(state, wave, scst) %>% 
  summarise("Electricity" = mean(elec_share),
            "Cooking Fuel" = mean(fuel_share)) %>% 
  gather(energy, share, -c(state, wave,scst)) %>% 
  ggplot(aes(scst, share, fill = energy)) +
  geom_col(position = "stack") +
  geom_text(aes(label = scales::percent(share, accuracy = 1)),
            position = "stack", vjust = 1) +
  facet_grid(wave ~ state) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_fill_brewer(palette = 5, type = "qual") +
  labs(y = "Mean share of monthly expenditure",
       x = NULL,
       fill = "Energy Use") +
  theme_bw() +
  theme(strip.background = element_rect(fill = NA, color = "black"),
        legend.position = "top")

ggsave(here("Manuscript","Figures","descr_energyburden.png"),
       device = "png", width = 8, height = 6*0.618, units = "in")

# Perception of stealing -------------------------------------------------------

pool %>% 
  transmute(state = case_when(
              state == 9 ~ "Uttar Pradesh",
              state == 10 ~ "Bihar",
              state == 19 ~ "West Bengal",
              state == 20 ~ "Jharkhand",
              state == 21 ~ "Odisha",
              state == 23 ~ "Madhya Pradesh"), 
            wave = ifelse(wave == 0, "2014", "2018"),
            others_steal = ifelse(is.na(others_steal), "No Answer", as.character(others_steal)),
            ) %>% 
  group_by(state, wave, others_steal) %>% 
  dplyr::count() %>% 
  group_by(state,wave) %>% 
  mutate(prop = round(n / sum(n),3)) %>% 
  ggplot(aes(wave, prop, fill = fct_rev(factor(others_steal, levels = c("No Answer", "Yes", "No"))))) +
  geom_col(position = "stack") +
  geom_text(aes(label = scales::percent(prop, accuracy = 1)),
            position = "stack", vjust = 1) +
  facet_grid(. ~ state) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_fill_brewer(palette = 5, type = "qual") +
  labs(y = "Proportion of respondents",
       x = NULL,
       fill = "Think co-villagers steal electricity") +
  theme_bw() +
  theme(strip.background = element_rect(fill = NA, color = "black"),
        legend.position = "top")

ggsave(here("Manuscript","Figures","descr_energythef.png"),
       device = "png", width = 8, height = 6*0.618, units = "in")
  